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I. INTRODUCTION 

Quantum tensor networks are a competitive tool to 
study strongly correlated quantum systems on a lattice. 
They originate from the density matrix renormalization 
group [1 - an algorithm to minimize the energy of a ma- 
trix product state (MPS) ansatz in one dimension (ID). 
In recent years the MPS was generalized to a 2D "ten- 
sor product state" better known as a projected entangled 
pair state (PEPS) [2J. Another type of tensor network 
is the multiscale entanglement renormalization ansatz 
(MERA) [3] that is, in some respects, a refined version 
of the real space renormalization group. Both PEPS and 
MERA can be applied to strongly correlated fermions 
in 2D [4 , because they do not suffer from the notorious 
fermionic sign problem. This makes them a powerfull tool 
to attack some of the hardest problems in strongly cor- 
related electronic systems, including the enigmatic high 
temperature superconductivity [5 j. Indeed, PEPS has al- 
ready provided first results for the ground state energy 
of the t — J model [6], that can compete with the best 
variational Monte-Carlo results [7]. 

In contrast to the ground state, thermal states have 
been explored mainly with the MPS in ID [8 , 9 , but they 
are more interesting in 2D, where they can undergo finite 
temperature phase transitions. In 2D thermal states were 
represented by tensor product states and contracted with 
the help of the higher-order singular value decomposition 
in Ref. [10]. A similar projected entangled-pair operator 
(PEPO) ansatz was proposed in Ref. [TTj . 

In this paper we follow a different route. In a way that 
can be easily generalized to 2D, the MPS can be extended 
to finite temperature by appending each lattice site with 
an ancilla [8 j. A thermal state is obtained by imaginary 
time evolution of a pure state in the enlarged Hilbert 
space, starting from infinite temperature. Unfortunately, 
in contrast to ID, where the time evolution of a MPS 
can be simulated accurately and efficiently, in 2D the 
time evolution of PEPS appears to be a hard problem. 
It requires accurate computation of a tensor environment 
that is often hard to approximate accurately and reliably. 
The aim of this work is to overcome this problem. 




FIG. 1. In A, graphic representation of the tensor A™ hl . In 
B, the amplitude with all bond indices connecting 

nearest-neighbor tensors contracted. The index contraction 
is represented by a line connecting two tensors. 

II. THERMAL STATES 

We consider spins on an infinite square lattice with a 
Hamiltonian H. Every spin has S states i = 0, S — 1 
and is accompanied by an ancilla with states a = 
0, ...,5 — 1. The enlarged Hilbert space is spanned 
by states risKs> a s)> where the product runs over lat- 
tice sites s. The state of spins at infinite temperature, 

p{(3 = 0) = n s (| 12i=o !j is obtained from 

a pure state in the enlarged Hilbert space, 

p(0) = T¥ ancillas |V>(0))(V>(0)| , (1) 

where 

w°)> = n(E^i^>) ( 2 ) 

is a product of maximally entangled states of every spin 
with its ancilla. The state p{f3) oc e~^ n at finite f3 is 
obtained from 

cc e-^w |V>(0)> = U{P) |V(0)> (3) 
after imaginary time evolution for time /3 with ^T~L. 

III. PEPS 

In the quantum Ising model with spin-^ that we con- 
sider in the rest of this paper, the translational invariance 
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is not broken and a unit cell encloses only one lattice site. 
Therefore, for an efficient simulation of the time evolu- 
tion we represent \ip(/3)) by a translationally invariant 
PEPS with the same tensor A^S^(/3) at every site. Here 
z,a = 0,...,5— 1 are the spin and ancilla indices respec- 
tively, 5 = 2, and r, 6, Z = 0, ...,D — 1 are the bond 
indices to contract the tensor with similar tensors at the 
nearest neighbor sites, see Fig. [I]A The ansatz is 
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Here the sum runs over all indices i s ,a s at all sites s. 
The amplitude is the tensor contraction in Fig. [1J>. 
The initial state Q can be represented by a tensor 
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D = 1 is the minimal bond dimension sufficient to repre- 
sent the initial state. 



IV. QUANTUM ISING MODEL IN 2D 

We proceed with 

n = - ^Z s Z s ,-h^X s = Uzz+Ux. (6) 

(s,s') s 

Here Z, X are Pauli matrices. The model has a ferro- 
magnetic phase with non-zero spontaneous magnetiza- 
tion (Z) for small h and large (3. At h = the critical 
point is f3 c = -ln(\/2 - l)/2 = 0.441, and at = 
the quantum critical point is h c = 3.04, see Ref. [13]. 



V. SUZUKI-TROTTER DECOMPOSITION 

We define U Z z(&P) = e~i nzzA(3 and U x (Af3) = 
e -2^xA/3 £ Qr interaction and the transverse field re- 
spectively. In the Suzuki- Trotter decomposition a small 
time step 

U(d/3) = U x (df3/2)U zz (df3)U x (df3/2) + (D(df3 3 ). (7) 
The action of U x {d0) on PEPS replaces A\ a rhl with 



oc A 



trbl 



3=0,1 



trbl 



(8) 



of the same bond dimension D. Here e = tanh d/J). 
However, the action of Uzz(df3) maps A to a new tensor 
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Here e = tanh(^d/3), indices St, s r , 85, 5^ G {0,1}, and 
s = St + 5 r + 55 + s/. This is an exact map, but B has 
the bond dimension 2D instead of the original D. 



VI. TENSOR RENORMALIZATION 

The bond dimension has to be truncated back to D 
in a way least distortive to the new PEPS \i/jb)- The 
general idea is to use an isometry W that maps from 2D 
to D dimensions: 
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see also Fig. |4p. The isometry should be the least de- 
structive to the norm squared (^b\^b)- The construction 
of the optimal isometry described in Figs. 2|3|4 is a vari- 
ant of the corner matrix renormalization [12] . It requires 
calculation of a tensor environment of B in the network 
representing (iPbI^b)- Unfortunately, this environment 
cannot be calculated exactly in an efficient way. This is 
why it is replaced by an effective environment, made of 
environmental tensors C, V, H, that should appear to the 
tensor B the same as the exact one as much as possi- 
ble. The environmental tensors are contracted with each 
other by indices of dimension M. Increasing M should 
make the effective environment more accurate. The over- 
all cost of renormalizing B back to the bond dimension D 
is polynomial in both D and M. It is dominated by the 
calculation of V in Fig. [3] that scales like M 3 D 4 when 
M > D A or M 2 D 8 otherwise. 

At the beginning of the evolution the environmental 
tensors C, V, H are initialized with random numbers and, 
in principle, they should be reinitialized after every time 
step. This, however, would not be the most efficient 
method for a smooth time evolution where both A and 
the environmental tensors change infinitesimally in an 
infinitesimal time step. Thus after every time step it 
might be more efficient to use the converged environ- 
mental tensors as the initial ones for the next step. This 
"recycling" would accelerate convergence in the next step 
because there would be very little to converge. However, 
we found that such "recycled" evolution is very fast in- 
deed but, especially near a phase transition, the tensors 
get trapped in lower dimensional subspaces, not making 
full use of the available dimension M. Results often ap- 
pear converged in increasing M while they are actually 
just trapped in an M e ff < M. This is not quite sur- 
prising, because even though the tensor A may evolve 
smoothly across a critical point, the environment does 
not need to be smooth, because it represents the rest of 
the infinite system at criticality. It is the environment 
that is critical, not A, even though the environment is 
made of an infinite number of smooth A's. To prevent 
the trapping, but at the same time not to slow down the 
algorithm too much, we add weak noise to the converged 
tensors before they are reused in the next time step. In 
practice, a noise at the level of < 1% of a typical ten- 
sor element was enough for the algorithm to produce the 
same results as if the tensors were reinitialized, but at a 
much faster rate. Since we want an accurate time evolu- 
tion, it is essential that the tensors do not get trapped in 




FIG. 2. In A, the contraction of two tensors B on the left 
hand side (LHS) gives the transfer matrix b on the right hand 
side (RHS), as seen from the top. In B, the contraction of 
transfer matrices on the LHS is the norm squared Tr p(/3) — 
| - This contraction cannot be done exactly. It is 
approximated by the contraction on the RHS with a corner 
matrix C and vertical/horizontal tensors V, H. Their (red) 
environmental indices have dimension M. The C, V, H are 
such that to the transfer matrix b in the center its environment 
on the RHS should appear the same as its exact environment 
on the LHS. Their construction is described in Fig. [3] 




FIG. 3. The tensors C,V,H are obtained by repeating a 
renormalization procedure until convergence. The procedure 
has two steps. In the first step the tensors C and H are 
contracted to form a matrix C.H. The M x M(2D) 2 matrix 
C.H is subject to singular value decomposition. It has M 
right singular vectors that define an isometry Z. The isometry 
is used to compress the right index of C.H back to dimension 
M giving a new corner matrix C' . The same Z renormalizes 
the contraction V.b giving a new V' . The second step is the 
same but with the roles of H and V interchanged. The two- 
step procedure is repeated until convergence measured by the 
figure of merit explained in Fig. [4] 

any single time step, because the errors can accumulate 
and derail the following evolution. 

Another technical issue concerns the construction of 
V in Fig. [3j In principle, all singular vectors Z of the 
corner matrix should be used in this contraction, even 
those corresponding to singular values equal to numeri- 
cal zero. The "zero vectors" do not make any difference 
when V' is contracted with C . However, we found the 
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FIG. 4. In A, for the isotropic tensors A the diagrams 
on the LHS would be two equivalent representations of the 
norm squared of the state, if not for the one uncontracted 
bond in the middle of each of them. For better numerical 
stability, we add these equivalent diagrams making a square 
matrix E of dimension {2D) 2 on the RHS. By construction, 
E is non-negative and its trace is equal to the norm squared. 
(Renormalization of a complex and non-symmetric E is de- 
scribed in Appendix pi) Tr E is the figure of merit refered to 
in the caption of Fig .^3] In B, each of the two indices of E can 
be represented by two indices of dimension 2D in such a way 
that the upper (lower) index corresponds to the top (bottom) 
layer of tensors B in Fig. [2|V. After the lower index is traced 
out we obtain a non-negative 2D x 2D matrix. Its D leading 
eigenvectors corresponding to the D largest eigenvalues define 
the isometry W . In C, the isometry renormalizes the new ten- 
sor B back to a tensor A' with the original bond dimension 
D thus completing the action of Uzz on PEPS. 



algorithm unstable unless we set the zero vectors in V 
to zero. These (numerically inaccurate) vectors do make 
a difference when V is contracted with a tensor other 
than C and this opens room for the observed instability. 



VII. ZERO TRANSVERSE FIELD 

In this classical limit the exact state = 
Uzz(P)\^(0)) can be obtained from the initial state by 
just one application of Uzz(P)- As in Eq. ([9|, this exact 
transformation doubles the bond dimension of the initial 
tensor Q to D = 2. Thus D = 2 (or D = S in general) 
is enough for an exact PEPS representation of any classi- 
cal state including the critical one. However, calculation 
of expectation values requires an approximate environ- 
ment build with the tensors C, V, H of a finite dimension 
M. The closer to criticality the bigger M is needed to 
represent the critical correlations in the environment. 

Figure [5] shows numerical simulations of the evolution 
by a product of small transformations Uzz(dP)- After 
each transformation the tensor is renormalized back to 
D = 2. The plots show excellent agreement with On- 
sager's solution, except in a narrow neighborhood of the 
critical point, but even there increasing M seems to con- 
verge the numerical solution towards the exact one. 
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FIG. 5. Numerical results versus exact Onsager's solution in 
the classical case of zero transverse field h = 0. Here (ZZ) is 
the ferromagnetic correlator between nearest neighbors, and 
(Z) is spontaneous magnetization. The upper panels show 
these quantities in a wide range of /3, and the lower panels 
zoom within 10 -3 of the critical /3 C . PEPS with D — 2 is an 
exact representation of a thermal state, but accurate evolution 
and calculation of expectation values near /3 C require large M. 
The lower panels show their convergence to the exact solution 
with increasing M. Here df3 = 10 _4 /3 C (upper panels) and 
d/3 = 10~ 6 /3 C (lower panels). 

The numerical results suggest that at the critical point 
a very large, if not infinite, M is needed for an accurate 
description of the long range critical correlations. Con- 
sequently, imaginary time evolution with a finite-M is 
bound to accumulate unrecoverable errors near the criti- 
cal point that will distort the following low temperature 
phase. To avoid this distortion, we suggest to add a tiny 
symmetry-breaking term to the Hamiltonian in order to 
smooth out the non-analyticity of the critical point and 
turn it into a smooth crossover. A PEPS with a finite- 
M can be evolved accurately across the crossover and 
into the low temperature phase. Event ualy the critical- 
ity can be recovered by turning the symmetry-breaking 
term down to zero and increasing M at the same time. 
This is what we do below in the quantum case of a finite 
transverse field. 

The spontaneous symmetry breaking in Fig. [5] may 
deserve a comment. In Eq. ^ the zero temperature fer- 
romagnetic state Uzz (00)1^(0)} is represented exactly by 
B™ Sr Sb Si oc (— iy s 5 za that does not break the symme- 
try. Its transfer matrix is bs^s^s^Si + 
jStjSrjSbjS^ where I± = (1,±1) T is a vector, j| is the 
5-th component of the vector I± : S t = s t + s t \ mod 2, 
and s's(s's) are the bond indices of the top(bottom) B in 
Fig. [2^. The 7±-part of b corresponds to (Z) = ±1, 
but the symmetry between these two parts is broken 
by the iterative construction of the environmental ten- 
sors. Indeed, for M = 1 the iterative procedure has 
two stable fixed points: C\\ — 1,^1,5,1 = H\s,i = I± 
breaking the symmetry to (Z) = ±1. The same is true 
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FIG. 6. The correlator Cr = (Z s Z a t) as a function of the 
separation R = \s — s'\ at different distances from the critical 
point measured by e = (f3 c — /3)//3 c . The log-log plot shows 
convergence to the exact function Cr ~ i? _1//4 (dashed line). 
Here D = 2, M = 32, and d{3 = 10" 6 /3 C . 

for M = 2 when the stable symmetry-breaking points 

are C Su s 2 = /f 1 ^ 2 , V Sl ,s„s 3 = H Sl ,s„s 3 = I^I^lt- 
By a suitable change of basis, each of these two fixed 
points can be represented by more compact tensors with 
M = 1. Thus the symmetry breaking reduces the re- 
quired M from 2 to 1. Once the symmetry is broken to 
a fixed point of the environment, the tensor renormal- 
ization in Fig. [4p also breaks the symmetry of the new 
renormalized tensor A'. This simple example explains 
the property of the algorithm observed in the ferromag- 
netic phase: the symmetric state is unstable but, once 
the symmetry is broken, the broken state is more accu- 
rate than the symmetric one for the same M. The broken 
state is simply less entangled than the symmetric one. 

VIII. FINITE TRANSVERSE FIELD 

For h > the Hamiltonian (|6| is quantum and a PEPS 
with a finite D is in general not an exact representation 
but an approximation to a thermal state. However, at 
finite temperature the fixed point of the renormalization 
group is a classical Hamiltonian whose critical thermal 
state can be represented by a PEPS exactly. This is 
why we expect that even a PEPS with the minimal non- 
trivial D = 2 (D = S in general) can in principle capture 
the universal critical properties of a quantum system at 
finite temperature, even though it may be not an accurate 
description of its short range quantum correlations. Just 
as in the classical case, it is mainly M and not D that 
limits the accuracy of PEPS at the critical point. 

In order to smooth out the finite-M imaginary time 
evolution across a critical point we added a small 
symmetry-breaking perturbation SH = —hz J2 S with 
a tiny longitudinal field hz- The perturbation is rounding 
the non-analyticity at the critical point making it possi- 
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FIG. 7. The magnetization (Z) as a function of f3 for the 
quantum Hamiltonian with h = \h c at different symmetry- 
breaking fields /iz — >• 0. For D = 6 (solid lines), D = 4 
(dashed lines), and D — 2 (dotted lines) the magnetization 
converges to a non- analytic critical curve when hz — »• 0. All 
plots are converged in M. The required M < 16 increases 
with decreasing /iz- An adaptive imaginary time-step d/3 was 
used with the shortest d/3 > 10~ 6 /3 C near the critical point. 

ble to evolve across the point with a finite M without 
accumulating unrecoverable errors. We expect that ac- 
curate evolution will require increasing M as hz is turned 
down to 0. 

Figure [7] shows numerical results for the magnetization 
(Z) at a relatively strong transverse field h = \h c . There 
is not much quantitative difference between the three sets 
of plots with D = 2,4,6. As expected, all three sets, 
even the minimal non-trivial D = 2, converge to a non- 
analytic critical curve when hz — > 0. Encouraged by the 
cross-section in Fig. [7| we also made a dense scan of the 
whole h — j3 phase diagram with the minimal D — 2. 
The ferromagnetic phase at low temperature and weak 
transverse field can be clearly read from the 3D plot in 
Fig. H 

IX. CONCLUSION 

A PEPS with ancillas can be efficiently evolved in 
imaginary time generating a PEPS representation of 
thermal states. In the case of a classical system, the 



bond dimension D equal to the number of states per 
site is enough for an exact representation of any ther- 
mal state. The evolution was simulated with the Suzuki- 
Trotter decomposition accurate to the second order in the 
time step. A variant of the corner matrix renormaliza- 
tion was used to obtain an accurate tensor environment. 
After every time step, the environmental tensors were 
perturbed by a weak noise to ensure that they make full 
use of their dimensionality. 

With some modification the algorithm can also evolve 
pure and thermal states in real time and, after intro- 
ducing fermionic swap gates, generate finite temperature 




FIG. 8. The magnetization (Z) as a function of the trans- 
verse magnetic field h and the inverse temperature /3. This 
3D plot clearly shows the ferromagnetic phase at low tem- 
perature and weak transverse field. Here D = 2, M = 24, 
h z = 10~ 10 , and d/3 = 10~ 2 /3 C . 



phase diagrams of strongly correlated fermions on a lat- 
tice [T4] . 
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FIG. 9. In A, the matrix E arises from the tensor net- 
work representing the norm-squared (^b\^b) after cutting 
the bond between the nearest-neighbor PEPS tensors B\ and 
B2. In B, the indices of E are renormalized by the projectors 
Pu and Py. In C, since the indices of E are also indices of 
the PEPS tensors B\ and B2, the PEPS tensors' indices are 
also renormalized by the projectors Pu and Py. The product 
PjjPv = YZ,p=i \ U <») ( u <*\ v p) ( V p\ inserted in the bond be- 
tween B\ and B2 results in a bond matrix S a p = {U a \Vp) on 
the bond B\ — B2. The singular value decomposition of the 
bond matrix S — u fi followed by absorbtion of the matrix 
u y/JX into the right index of B\ and the matrix y/JX into 
the left index of B2 completes renormalization of the bond 
Bi — B2. After renormalization of all bonds a new PEPS is 
obtained with new tensors A' . 



Appendix A: renormalization of general complex 
and non-symmetric PEPS tensors 

The matrix in Fig. [4^, we will call it E here, is 
used for the renormalization of the PEPS tensors from 
B to A! . Its trace is the norm-squared of the PEPS: 
TtE = (^b\^b)- The matrix itself is the norm-squared 
in Fig. [2j3, but with one of the bonds connecting pairs 
of nearest-neighbor tensors B cut open. For the isotropic 
tensors considered in this paper E is by construction real, 
symmetric, and positive semi-definite. However, in gen- 
eral a PEPS may be neither isotropic, nor translation- 
ally invariant, nor even real, and the matrix Eij is just a 
complex matrix. We show such a general matrix E for a 
bond between inequivalent PEPS tensors B\ and B 2 in 
Fig. |9]a. We want to renormalize the indices of this E, 
because they are also the indices of the PEPS tensors B\ 
and B2 that have to be renormalized back to the original 
bond dimension D. 

The renormalization procedure begins by a singular 



value decomposition 

2D 

e = Y.\ U ^ X ^\ - ( A1 ) 

Here A's are the singular values in decreasing order 
Ai > A2 > and \U a ) (\V a )) are corresponding left 
(right) singular vectors. We define projectors Pjj = 

Ea=i Pa)(U a \ and P v = £f =1 \V a )(V a \. A renormal- 
ized matrix is 

D 

& = P v EPy = ^2\U a )X a (V a \ , (A2) 

where we truncate the sum ( |A1[ ) from 2D to D, see Fig. 
This truncation minimizes the difference between E 
and the renormalized E' measured by the error 

2D 

Tr(E — E r )\E — E r ) = X l ' ( A3 ) 

a=D+l 

As mentioned above and shown in Fig. [9]A., the indices of 
E that are renormalized by the projectors Pjj and Py are 
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at the same time the bond indices of the nearest-neighbor 
PEPS tensors B 1 and B 2 , see Fig. [9p. 

The expression Tr£" is the norm-squared of the orig- 
inal PEPS, but with an additional "bond matrix" S in- 
serted in the bond connecting the renormalized nearest- 
neighbor PEPS tensors: 

S af3 = (U a \V p ) . (A4) 

Here a, /3 = 1, ...,£). In order to go back to the original 
PEPS structure, without any bond matrices, we want to 
absorb the bond matrix into the PEPS tensors connected 
by the bond. To this end, we make one more Schmidt 
decomposition 



Here u : v are unitary D x D matrices and /i is a diag- 
onal matrix of singular values. In Fig. [9p the matrix 
u ^fjl is absorbed to the left PEPS tensor B\, and the 
matrix ^JJi to the right tensor £? 2 , thus completing 
the renormalization procedure. Alternatively, the whole 
bond matrix S can be simply absorbed to any of the two 
PEPS tensors connected by the bond. We use the more 
symmetric version to simulate evolution of PEPS in real 
time [H]. 



S = U /! V 



(A5) 



